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Abstract 

We present the problematic of controlling the discreteness effects in cosmological N-body sim- 
ulations. We describe a perturbative treatment which gives an approximation describing the 
evolution under self-gravity of a lattice perturbed from its equilibrium, which allows to trace 
the evolution of the fully discrete distribution until the time when particles approach one an- 
other ("shell-crossing"). Perturbed lattices are typical initial conditions for cosmological N-body 
simulations and thus we can describe precisely the early time evolution of these simulations. A 
quantitative comparison with fluid Lagrangian theory permits to study discreteness effects in 
the linear regime of the simulations. We show finally some work in progress about quantifying 
discreteness effects in the non-perturbative (highly non-linear) regime of cosmological N-body 
simulations by evolving different discretizations of the same continuous density field. 
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1. Introduction 

In this proceedings we present some specific work about the material presented in the 
proceedings of Michael Joyce ^ , where a detailed introduction can be found. In what 
follows we remind very briefly the motivations of our work. 

At the scales of interest in cosmology, the evolution of the dark mater can be very well 
described as a fluid, which obeys the Vlasov equation coupled with the Poisson equation 
(e.g. (1)): 
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where /(r,v,t) is the one-particle phase space density in physical coordinates (i.e. it 
gives the probability density to find a particle at r with velocity v at time t), p{r,t) 
the local density of dark matter and G the gravitational constant. The solution of the 
Vlasov equation is generally "estimated" performing N-body simulations, using a number 
of particles much smaller than the number of dark matter particles (it can be seen as a 
discretization of the system). This method leads inevitably to discreteness effects, which 
are still not well understood. 

Up to now the primary approach to the study of discreteness in N-body simulations 
has been through numerical studies of convergence (see e.g. (2; 3)), i.e., one changes the 
number of particles in a simulation and studies the stability of the measured quantities. 
Where results seem fairly stable, they are assumed to have converged to the continuum 
limit. While this is a coherent approach, it is far from conclusive as, beyond the range of 
perturbation theory, we have no theoretical "benchmarks" to compare with. Nor is there 
any systematic theory of discreteness effects, e.g., we have no theoretical knowledge of the 
N dependence of the convergence. It can be found in the literature several studies using 
this approach of the most obvious effect of discreteness — and certainly the one most 
emphasized in the literature — which is the two body coUisionality: pairs of particles can 
have strong interactions with one another, which is an effect absent in the coUisionless 
limit. For analysis and discussion of these effects see, e.g., (4; 5; 6; 7; 8). 

In this proceedings we present a work which permits to understand systematically the 
discreteness in the linear regime of gravitational clustering. We do so by developing a 
perturbative solution to the fully discrete cosmological N-body problem. This essentially 
analytic solution can be compared to the analogous fluid (N —>■ oo) theory, of which 
analytic solution is well-known. We will remind first the standard perturbative fluid theory 
(FLT), we will construct the particle linear theory (PLT) and by comparing both we will 
quantify the discreteness effects in the linear regime. We will also present a test which 
permits to quantify the discreteness effects in the non-perturbative regime of gravitational 
clustering. The material presented in these proceedings is based on work which can be 
foimd in (9; 10; 11). 

2. Fluid theory 

It is convenient, to simplify the problem Vlasov equation, to construct a set of fluid 
equations, which give a less detailed (but sufficiently accurate) description of the system. 
We define the mass density and the mean fluid velocity from the velocity moments of 



where m is the mass of the dark matter particles. Inserting Eqs. (2) in the Vlasov-Poisson 
equation (1), we obtain the continuity equation and the Euler equation in the variables 
p(r, t) and v(r, t). These are the fluid equations in the Eulerian formalism. In cosmology, 
it is common to use the Lagrangian formalism which is known to give better results 



/(r,v,t): 





2 



in a pcrturbativc approach (see e.g. (12)). The idea of the Lagrangian formulation is to 
calculate the trajectories of infinitesimal fluid elements. The velocity is given by the 
velocity of these fluid elements and the density varies according with the convergence or 
the divergence of the fluid elements to each point. 

We define the Lagrangian coordinate q as the position of the fluid element at the 
initial time ^ . Because we work in an expanding universe, it is natural to write the 
physical position of a fluid element r in function of its Lagrangian coordinate q as: 

r(i) = a(t)(q + u(q,t)), (3) 

where a{t) is the scale factor (which is given by the particular cosmological model consid- 
ered) and u(q, t) is a "displacement field" of the fluid elements. The displacement field 
u(q, t) in Eq. (3) can be understood as a perturbation of the homogeneous model which 
expands with the scale factor a(t). Using the coordinate transformation (3) we obtain 
a set of fluid equations for the displacement field u(q, t). Linearizing them about the 
displacements u(q, t) and neglecting pressure corrections (which is a very good approxi- 
mation for sufficiently large scales), one obtains a simple system of equations (sec e.g. (13) 
for details). To look for a solution of these equations, it is convenient to divide the dis- 
placements into a curl-free part, U[| and a divergence-free part, uj^, i.e., u = uy -|-u^ (i.e. 
V X U|| = and V ■ = 0). Then we obtain the set of equations (choosing appropriate 
boundary conditions) (13): 

Ui+2-Ui=0, iiii +2-U|| -h3-U|| = 0. (4) 

Given initial conditions, it is simple to compute the solution for the displacement field 
with time. The linear approximation (and more generally a perturbative approach of the 
Lagrangian fluid equations) breaks down when the volume of a fluid element becomes 
zero. This is called sheJl crossing. 

3. Linear perturbative theory 

In this section we will present a linearized theory of the N-body problem (PLT), which 
is the discrete analogue of the linear fluid theory presented in the previous section. The 
starting point is the full evolution of the N-body system in comoving coordinates: 

a _ 1 ^ Gm(x, - Xj) 

where x^ is the position of the i-tli particle and a{t) the scale factor. 

A standard way to generate initial conditions for N-body simulations consists in per- 
turbing a lattice (see e.g. (14)). It is therefore natural to build a perturbative theory 
where the perturbed variable is the displacement of each particle about the lattice, which 
is an equilibrium position. We will therefore have an accurate description of the clustering 
when the displacements (or, in fact, the relative displacements) are smaller than the inter- 
particle distance. When the relative displacements become larger than the inter-particle 
distance, the approximation breaks down, which is equivalent to the "shell-crossing" in 
fluid theory. 

^ It can be viewed just as a "label" of the particle. 
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Figure 1. Spectrum of eigenvalues for simple cubic lattice with 16^ particles. The lines correspond to 
chosen directions in k space. 

Wc define the displacement u(R, t) about the lattice position R as Xi(t) = R+u(R, t), 
where R is the lattice site of particle i. Using this notation we can linearize the full 
evolution (5) as 

u(R,t) + 2-u(R,i) = -XT^VrR,--R')vL{Ii',t), (6) 

R' 

where I'(R) is called the dynamical matrix, which is ^'^^(R = 0) = — X^R^i^o ■^a"' (-f^) 
and X'p;^(R ^ 0) = Gm — 3 j . Taking periodic boundary conditions (as in the 
N-body simulations), it is possible to diagonalize Eq. (6) simply by taking its Fourier 
transform, defined as u(k, i) = u(R, t)e~'*''^ for the displacement field, and in an 
analogue manner for the dynamical matrix I'(k). Using Eqs. (6) and the definition of 
the Fourier transform, we obtain for each k: 

Q(k, t) + 2-Q(k, t) = -^P(k)u(k, t). (7) 

a a-^ 

It is possible to solve Eq. (7) by diagonalizing the 3x3 matrices 2'(k). For each k, this 
determines three orthonormal eigenvectors e„ (k) with three associated real eigenvalues 
Ci;^(k) {n = 1,2,3), satisfying the eigenvalue equation 2?(k)e„(k) = ti;^(k)e„(k). Once 
a;^(k) is calculated, the evolution of the displacement of the particles is straightforward 
to compute. 

3.1. Spectrum of eigenvalues o/X'(k) 

To understand the dynamics of the PLT is important to study the spectrum of eigen- 
values of the matrix 2?(k). We will describe in what follows this spectrum for a simple 
cubic (hereafter sc) lattice, which is very widely used in N-body simulations of structure 
formation in cosmology. 

In Fig. 1 wc plot the spectrum of a sc lattice, for N = IG'^ particles. We show the 
normalized eigenvalues e,i(k) = (k)/47rG(Oo as a function of the modulus of the k 

^ The eigenvalues are real and the eigenvectors orthonormal because 'D(k) is a real and symmetric 
matrix. 
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Figure 2. Amplification function A'^(\i,t) divided by the fluid amplification factor at a = 5, for a sc 
lattice, in function of the modulus of the wave-vectors normalized to the Nyquist frequency. 

vectors, normalized to the Nyquist frequency ~ tt/£ {£ is the intcr-particlc distance). 
With this normahzation the spectrum remains substantially the same as we increase the 
number of particles: the only change is that the eigenvalues become denser in the plot, 
filling out the approximate functional behaviors with more points. We see how for each 
vector k there are three eigenvalues uj'^{'k), n = 1,2,3. Each family of eigenvalues (i.e. 
with same n) defines a surface, corresponding to the three branches of the frequency- 
wavevector dispersion relation. Sections of these surfaces are plotted for some chosen 
directions of the vector k in Fig. 1. 

The fluid limit of this system is given by taking the k ^ limit keeping the inter- 
particle distance i constant: a plane wave fluctuation e^^'^ with k <C 1/^ has a variation 
scale much larger than the inter-particle distance, and therefore does not "see" the parti- 
cles. It is simple to show (e.g. (10; 15)) that in this limit the eigenvectors and eigenvalues 
arc: (i) one longitudinal eigenvector polarized parallel to k with normalized eigenvalue 
ei(k — > 0) = 1 and (ii) two transverse eigenvectors polarized in the plane transverse 
to k with normalized eigenvalues £2,3 (k ^ 0) = 0. It is interesting to note that in this 
limit, using Eq. (6), we obtain exactly the two fluid equations (4), particularized to an 
Einstcin-dc Sitter (EdS) universe and with the same boundary conditions. This means 
that the N-body system is a specific discretization — at least in the linear regime — of 
the underlying fluid theory. 

3.2. Discreteness effects 

It is now straightforward to quantify precisely the discreteness effects — in the linear 
regime — by comparing the FLT explained in section 2 and the PLT derived in section 
3. We will focus here on the discreteness effects in the power spectrum (hereafter PS, 
for its definition see e.g. (1; 16)) but an analogous analysis can be performed for any 
statistical quantity or, for example, the trajectories of a single particle (see (10)). 

It is possible to show, considering an initial power spectrum Po(fc) at t = tpj that we 
have, as a very good approximation (for not too large displacements compared to the 
lattice spacing), P(k, t) ~ j4p(k, t)Po{k), where j4p(k, t) is a function which depends on 
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Figure 3. Evolution of the PS and Ci(k) for, from up to bottom and left to right, a = 1, a = 2^ and 
a = 2^. We have also included the evolution of the PS predicted by linear theory ("LT"). 



the particular lattice and the cosmological model considered. One of the most important 
result obtained in this work is that discreteness effects increase with time. In the literature 
the time was not considered as a parameter in which depends the discrctness effects, but 
only "static" parameters as the number density of particles. 

It is possible to check that, in the fluid limit, we recover the well-known amplification 
of the PS in the fluid theory, i.e., limk_,o Ap{k, t) ~ {t/t^f'^^ = a. We plot the function 
j4p(k, t), normalized to the fluid amplification, in Fig. 2. We have chosen a value of 
a = 5 for the scale factor. This is a typical scale factor at which shell crossing occurs in 
cosmological simulations. Notice the similarity of this figure with the optical branch in 
Fig. 1: the evolution "deforms" the spectrum of eigenvalues. Note how the eigenvalues 
with e > 1 give rise to ylp(k, t) > a for these modes, i.e., there are modes which grow 
faster than the fluid limit. In the figure, we have classified the modes as a function of 
the angle subtended by their wave vector k with the lattice axis that form the minimal 
angle with it. We see that there is a strong dependence of the value of the eigenvalue on 
this angle: the closer k is to parallel to one of the axes, the larger is the eigenvalue of the 
mode, on average. This is a manifestation of the breaking of isotropy introduced by the 
N-body discretization on the lattice. Even if there are some modes that grow faster than 
the fiuid, averaging over bins with similar |k| the resultant growth is slower — because 
we consider sufficiently early times — than the fluid limit. Note that this averaging is 
generally performed when computing statistical properties of the particle distribution, 
such as the PS, for example. 



4. Discreteness effects in the non-perturbative regime 

It is much more difficult to quantify the discreteness effects in the non perturbative 
regime because there is no known solution of the fluid theory. In this section we present 
briefly some work in progress which goal is to quantify the discreteness effects in this 
regime. We run a set of simulations with exactly the same parameters (time-step, force 
accuracy, etc.) for an initial PS P{k, 0) oc (which corresponds to the kind of density 
fluctuations of the standard cosmological model) in an expanding EdS universe, taking 
different perturbed lattices to generate the initial conditions: body cubic centered (bee). 
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face cubic centered (fee), sc and "glass" . It can be considered as different discretizations 
of the same continuous density field. 

In Fig. 3 we show the evolution of the PS for times corresponding to a = 1, a — 2^^, 
a = 2^. In the largest figures arc shown the PS and in the smallest ones the relative 
residuals Ci{k) between the PS of a particular lattice i {i = bee, fee, glass or sc) and the 
average of all of them: 

1 mPdk) - V'" . PJk) 

At the initial time, the PS up to the Nyquist frequency is approximately the same for all 
the lattices. However, at scales larger than the Nyquist frequency, the PS corresponding to 
the different lattices is very different due to the great differences between the distributions 
at small scales. Time evolution (a = 2"^) erases almost all the differences at small scales. 
At all scales, the evolution of the PS of the different initial distributions is very similar. 
The residuals show that i) the discreteness effects have been globally reduced, ii) the 
differences can be up to 10% and Hi) the scale at which these differences appear is much 
larger than the initial discreteness scale kN'- there is a propagation of discreteness from 
small scales to large scales. In the figure we have shown also the evolution of the initial 
PS using fluid linear theory. At a = 2^ there are no linear modes left in the box (it can be 
seen by the fact that linear fluid theory does not reproduce the evolution of any smallest 
Fourier mode). The differences between the PS of the different initial lattices are larger 
than in the previous time slice a = 2"^ (in which there were still some linear modes) . This 
is related to the fact that all the system is highly non linear. Finally, it is interesting to 
note that, as in the linear regime, the discreteness effects increase with time. 
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